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Abstract 

A linear dispersive mechanism leading to a burst in the Loo norm of 
the error in numerical simulation of polychromatic solutions is identified. 
This local error pile-up corresponds to the existence of spurious caustics, 
which are allowed by the dispersive nature of the numerical error. From 
the mathematical point of view, spurious caustics are related to extrema 
of the numerical group velocity and are physically associated to inter- 
actions between rays defined by the characteristic lines of the discrete 
system. This paper extends our previous work about classical schemes to 
dispersion-relation preserving schemes. 
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1 Introduction 

The analysis and the control of numerical error in discretized propagation-type 
equations is of major importance for both theoretical analysis and practical ap- 
plications. A huge amount of works has been devoted to the analysis of the 
numerical errors, its dynamics and its influence on the computed solution (the 
reader is referred to classical books, among which [SJ [T31 [51 [S| ) ■ The emergence 
of Dispersion-Relation-Preserving (DRP) schemes [3]), which have the same dis- 
persion relation as the original partial difference equations, enables one to have 
very accurate high order finite difference schemes. 

The two sources of numerical error are the dispersive and dissipative properties 
of the numerical scheme, which are very often investigated in unbounded or 
periodic domains thanks to a spectral analysis. 
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It appears that existing works are mostly devoted to linear, one-dimensional 
numerical models, such as the linear advection equation 



du du 







(1) 



dt dx 



where c is a constant uniform advection velocity. 

The two sources of numerical error are the dispersive and dissipative proper- 
ties of the numerical scheme, which are very often investigated in unbounded 
or periodic domains thanks to a spectral analysis. Following this approach, a 
monochromatic wave is used to measure the accuracy of the scheme. Such a tool 
is very powerful and provides the user with a deep insight into the discretization 
errors. But some results coming from practical numerical experiments still re- 
main unexplained, despite the linear character of the discrete numerical model. 
As an example, let us note the sudden growth of the numerical error for long 
range propagation reported by Zingg [15] for a large set of numerical schemes, 
including optimized numerical schemes. 

The usual modal analysis is almost always applied to monochromatic reference 
solutions, with the purpose of analyzing the error committed on both their am- 
plitude and their phase, leading to classical plots of the relative error as the 
function of the Courant number and/or the number of grid points per wave- 
length. Therefore, dispersive phenomena associated to polychromatic solutions 
are usually not taken into account. 

The present paper deals with the analysis of linear dispersive mechanism which 
results in local error focusing, i.e. to a sudden local error burst in the L^o norm 
for polychromatic solutions. This phenomena is reminiscent of the physical one 
referred to as the caustic phenomenon in linear dispersive physical models [2] , 
and will be referred to as the spurious caustic phenomenon hereafter. It extends 
our previous work [2] to DRP schemes. The present analysis is restricted to 
interior stencil, and the influence of boundary conditions will not be considered. 
The paper is organized as follows. Main elements of caustic theory of interest 
for the present analysis are recalled in section [21 DRP schemes are presented in 
section [31 Their caustical analysis is exposed in section [31 A numerical example 
is presented in section [5l 

2 Caustics 

The solution of Eq. ^ is taken under the form: 



where w = + i 77^^ is the complex phase, and k the real wave number. For 
dispersive waves, it is recalled that the group velocity Vg{k) is defined as 



(2) 




(3) 
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A caustic is defined as a focusing of different rays in a single location. The 
equivalent condition is that the group velocity exhibits an extremum, i.e. there 
exists at least one wave number kc such that 




(4) 



The corresponding physical interpretation is that wave packets with character- 
istic wave numbers close to fc^ will pile-up after a finite time and will remain 
superimposed for a long time, resulting in the existence a region of high energy 
followed by a region with very low fluctuation level. 

The linear continuous model Eq. ([1]) is not dispersive if the convection velocity 
c is uniform, and therefore the exact solution does not exhibits caustics since the 
group velocity does not depend on k. The discrete solution associated with a 
given numerical scheme will admit spurious caustics, and therefore spurious local 
energy pile-up and local sudden growth of the error, if the discrete dispersion 
relation is such that the condition ([4]) is satisfied. For a uniform scale-dependent 
convection velocity, such spurious caustics can exist in polychromatic solutions 
only, since they are associated to the superposition of wave packets with different 
characteristic wave numbers. 



The general dispersion relation associated with the discrete scheme enables us 
to obtain the corresponding group velocity, given by: 



The numerical solution will therefore admits spurious caustics if 



The corresponding values of if and k will be respectively denoted (pc and fee- 
Spurious caustics are associated with characteristic lines given by 



Set: 




(5) 






(8) 



where 



(9) 



3 DRP schemes 



The Burgers equation: 



Ut + CUUx — 



fJ'Uxx = 0, 



(10) 
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c, fj, being real constants, plays a crucial role in the history of wave equations. 
It was named after its use by Burgers [Tj for studying turbulence in 1939. 

i, n denoting natural integers, a linear finite difference scheme for this equation 
can be written under the form: 

J2ai^ur = (11) 

where: 

ur = u{lh,mT) (12) 

I S {j — 1, i, « + 1}, m G {n — 1, n, n + 1}, j = 0, rix, n = 0, nt- The 
aim are real coefficients, which depend on the mesh size h, and the time step r. 
The Courant-Friedrichs-Lewy number (cfl) is defined as a ~ ct / h . 
A numerical scheme is specified by selecting appropriate values of the coeffi- 
cients aim- Then, depending on them, one can obtain optimum schemes, for 
which the error will be minimal. 



m being a strictly positive integer, the first derivative ^ is approximated at 
the Z*'' node of the spatial mesh by: 

(13) 

k— — in 

Following the method exposed by C. Tam and J. Webb in the coefficients 
gammak are determined requiring the Fourier Transform of the finite difference 
scheme to be a close approximation of the partial derivative ( §f 
is a special case of: 

E Ikuix + kh) (14) 

k— — m 

where a; is a continuous variable, and can be recovered setting x — Ih. 
Denote by w the phase. Applying the Fourier transform, referred to by ^ , to 
both sides of pi]) , yields: 

m 

juju- -fke^'"^^u (15) 

j denoting the complex square root of —1. 

Comparing the two sides of enables us to identify the wavenumber A of the 

rn 

finite difference scheme (|13p and the quantity i 7^ e '\ the wavenumber 

k— — m 

of the finite difference scheme ((T^ is thus: 
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A 



fc=- 



j k (jj h 



(16) 



To ensure that the Fourier transform of the finite diff'erence scheme is a good 
approximation of the partial derivative over the range of waves with wave- 

length longer than 4 h, the a priori unknowns coefficients 7fe must be choosen 
so as to minimize the integrated error: 

£ = J\\Xh~Xh\'^d{Xh) 



C - 7fc sm(fcC) 

fc — — 771 



= C- E 7fcSin(fcC) 

The conditions that £ is a minimum are: 



^ 7fc cos(fcC) 

; — — m 

^ 7fe cos(fcC) 



(17) 



(18) 



1. e.: 



I -Csin(iC) + 51 '^''^ cos((fc-i)C) WC = 
Changing i into — z, and k into — in the snmmation yields: 



(19) 



1. 6.: 



Thus: 



Csin(iC)+ 7_fe cos((-fc + i)C) WC = 



k— — m 



Csin(iC)+ 7-fc cos((fc-i)C) MC = 



(20) 



(21) 



Y ^^-k + 7fc} COS ( (fc - i) C) dC = 



k— — r. 



which yields: 



- {7_, + 7, 



a+ E { 

k—~-7n 



7-fc + 7fc 
k — i 



sin((fc-z)-j =0 



(22) 



(23) 
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which can be considered as a linear system of 2 m + 1 equations, the unknowns 
of which are the 7_j +7i , i = — m, . . . , m. The determinant of this system is not 
equal to zero, while it is the case of its second member: the Cramer formulae 
give then, for i = — m, . . . , to: 



or: 



For i = 0, one of course obtains: 



All this ensures: 



7_j + 7j = 



70 = 



(24) 
(25) 
(26) 



(27) 



k— — r, 



The values of the -fk coefficients are obtained by substituting relations ([25|l into 
(HI: 



(28) 



k— — r, 



TO being a strictly positive integer, a 2to + 1-points DRP scheme ([3]) is thus 
given by: 



(29) 



where the 7^, /c G {— m, to} arc the coefficients of the considered scheme, and 
satisfy the relations (p5|) . 



4 General study of DRP schemes 

The dispersion relation related to a general £>i?P-scheme is given by: 

m. 



h 



from which it comes that 



1 = 



(30) 



j T = B r — In 



E 



Ik e 



ik ifi^ 



1 + 



V 



(31) 



The group velocity can be expressed as 



6 



i k 



Va 



, m 



(32) 



7fc e 



i k ifi 



from which it comes that 



dip 



cr I c + cr > 7fe e 
fc— — m 



Through identification of the real and imaginary part of ([33]), we obtain: 



(33) 



a 7fc 7j_^^ {/c"^ sin [ {k + — fc/ cos [ {k + = — c fc^ 7^ sm{k(p) (34) 

k,l= — m k— — m 

and 

m m 

CT 7fc 7; {— COS [ (A; + /) y:] — A; / sin [ (A; + /) i^]} = c A;^ 7^ cos(A; t^) (35) 

k,l — — m k=~-m 

Due to (P5|) . ([35]) and ((37)) respectively become: 

m m 

cr 7fc 7^ |A;^ sin [{k -\- l)ip] — kl cos [ (A: + /)f^ ] } = —2 c A;"^ 7^ sin(A: (p) (36) 

k,l — — m k — l 

and 

m 

o" E 7fc 7; {" COS [(fc + OV'] - fc' sin [(fc + '/']} = (37) 

A; . / — — 7?i 

Denote by Tj, j e N* the Chebyshev polynomial of the first kind, and by Uj, 
j £ IN* the Chebyshev polynomial of the second kind: 



[f] 

cos(jx)=T,(cos(.t)) = ^^(-1) 



{n-k- 1)! 
k\{n-2k)\ 

sin(j x) ~ sin(x) Uj (cos(x)) 



(2 cos(.t)) 



n-2k 



k=0 



where: 



C/,(cos(x)) = ^(-l) 



k=0 



kl {n-2k)\ 



(2 cos(x)) 



(38) 
(39) 

(40) 



[^] denotes the integer part of J. 
Equations ((36)) . (jST]) can thus be written as: 



X] 7fe 7; {fc^sin(vp) i7fe+; (cos((p)) - fc«Tfe+; (cos((/3))} = -c ^ fc^ 7fe sin((/3) C/j; (cos((/3)) 

k,l — — rn h = — m 

(41) 

and 

o- 5Z (^sCi^")) + k I sin{^p) Uk+l (cos(i/p))} = (42) 

Using the relation: 



sin((p) = -cos2(^) (43) 
equations (|4T|) . ((42|) can be written as: 

h (cos(^)) = (44) 

and 

h {cos{^)) = (45) 

where, for all ^ G M: 

m m 

h{e) = a J2 -yk-yi {kWi-e^ Uk+l {9)- k I n+i (6)] + c ik Vi - (e) 

k,l— — m k — ~m 

(46) 

i.e.: 

m m 

/i(0) = (T ^ 7fc7; C/fc+ae)-fc«Tfe+ae)} +2c ^fcSfe Vl-92C/fc(e) 

k,l = — m, k = l 

(47) 

and 

/2(e) = ^ 7fc 7i {Tfc+i (9) + fcZ Vl - 02 C7fc+, (0)} (48) 

Due to: 

Tj (1) = 1 V j e N* (49) 

it is worth noting that: 

m 

/i(l) = -a ^ -tklikl (50) 

A: , / — — m 

and 

m 

/2(1)- 5] 7fc7i (51) 

J — — m 
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The knowledge of the scheme coefficients 7^, k E {— to,to}, enables one to 
study their variations, and to determine wether the equations (|44|) . (|45l) admit 
a solution. One can thus know wether 



the schema has spurious caustics. 



dip 



admits real roots, i. e. wether 



5 Numerical application: the 3-points DRP scheme 

The 3-points DRP scheme is given by: 

71 = 0.63662 



We thus have: 



and 



/i(l) = -2^{72-72}=0 



For the 3-points DRP scheme, the dispersion relation is: 



which leads to: 



It yields: 



+ e 



ir 



(-0.63662 + 0.63662 



^ — Arctan 

T 



e ^'^ 0.63662 cr(e2*¥' _ 1) 
(1 + 0.63662 cr) sm{(p) 



1 + (0.63662(7- 1) cos{ip) 



The derivative of the group velocity Vg vanishes for ip ^ 0, ip 



±7T 



(52) 
(53) 

(54) 

(55) 

(56) 

(57) 

and 



V? = 0.950935. 

The 3-points DRP scheme admits thus spurious caustics. 

Wc now illustrate the caustic phenomenon considering the two following sinu- 
soidal wave packets: 



Ml 



U2 



e (^-^o-ctf Cos [ fci (x - 4 - ct) ] 



(58) 



(59) 



where a > 0. The two wave packets are initially centered at Xq and Xq, re- 
spectively. The group velocity of the two wave packets arc Vi ~ Vg{ki) and 
V2 = Vg{k2), respectively, where the function Vg{x) is associated to the numer- 
ical scheme used to solve Eq. ([T]). 

If the solution obeys the linear advcction law given by Eq. ((T|) , the initial field is 
passively advected at speed c, while, if the advection speed is scale-dependent (as 
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in numerical solutions), the two packets will travel at different speeds, leading 
to the rise of discrepancies with the constant-speed solution. Another dispersive 
error is the shape-deformation phenomenon: due to numerical errors, the exact 
shape of the wave packets will not be exactly preserved. This secondary effect 
will not be considered below, since it is not related to the existence of spurious 
caustics. It is emphasized here that the occurance of spurious caustics originates 
in the differential error in the group velocity, not in the fact that shapes of the 
envelope of the wave packets are not preserved. The issue of deriving shape- 
preserving schemes for passive scalar advection has been adressed by several 
authors (e.g. [SUZ])- 

The spurious caustics will appear if the two wave packets happen to get su- 
perimposed. During the cross-over, the L^o norm of the error (defined as the 
difference between the constant-speed solution and the dispersive one) will ex- 
hibit a maximum. The characteristic life time of the caustic, t* , depends directly 
on the difference between the advection speeds of the two wave packets and the 
wave packet widths. Denoting li and I2 the characteristic length of the two wave 
packets, the time during which they will be (at least partially) superimposed 
can be estimated as 

It is seen that, since caustics are defined as solutions for which dVg/dk = 0, t* 
will be large if \ki — ^2! ^ 1- Noting fci ~ kc + 6k and k2 ~ k^ — Sk, one obtains 

m (61) 

2{Sk)^ ^(fc,) 

leading to t* oc (5fc)~^. 

Neglecting shape-deformation effects and assuming that the numerical scheme 
is non-dissipative, the numerical error E is given by: 



E = I e-"(^-^«-'=*)'Cos[fci {x-xl~ ct) ] - e-"(^-^o-*^i)'Cos[fci(a; - - tVi)] 

+ e-"(^-^«-"*)'Cos[fc2(a; - 4 - ct)] ~ (.-<^{^-A-iy^f Cos [fca {x - xl - t F2) 

(62) 

A simple analysis show that 



lim Loo{E{t)) = Loo{ui{t = 0)), maxL^{E{t)) = 2Loo{ui{t = 0)) (63) 

t — *+oc t 

The time history of the Loo norm of E for the 3-points DRP scheme scheme, is 
displayed in Fig. [TJ showing the occurance of the caustic and the sudden growth 
of the Loo error norm. 

Figure [5] displays the isovalues of the residual kinetic energy for 3-points DRP 
scheme, for cfl = 0.9. Minima are in black, maxima in white. In each case, the 
caustic corresponds to the white domain, where the residual kinetic energy is 
maximal. 
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Figure 1: Time history (Loo norm) of the numerical error for the two- wave 
packet problem (shape deformation and dissipative errors are neglected to em- 
phasize the linear focusing phenomenon). Numerical parameters are a = 0.0005, 
h = 0.01, Vi = —2.68381, V2 = —2.51381, corresponding to the properties of 
the 3-points DRP scheme, for a = 0.9. 




Figure 2: Isovalues of the residual kinetic energy for the 3-points DRP scheme, for 
cfl = 0.9. 

6 Concluding remarks 

In the above, we have set a general method that enables one to determine wether 
a DRP scheme admits or not spurious caustics. 

The existence of spurious numerical caustics in linear advcction DRP schemes 
has been proved. This linear dispersive phenomenon gives rise to a sudden 
growth of the Loo norm of the error, which corresponds to a local focusing of 
the numerical error in both space and time. In the present analysis, spurious 
caustics have been shown to occur in polychromatic solutions. The energy of 
the caustic phenomenon depends on the number of spectral modes that will get 
superimposed at the same time. As a consequence, the spurious error pile-up will 
be more pronounced in simulations with very small wave-number increments. It 
has been shown that a popular existing scheme, as the 3-points L'iJL'-scheme, 
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allows the existence of spurious caustics. 
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